knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)

TIDY DATA:

Cleaning up explan

First chunk of code subsets to only using 8 blocks per survey location, 4 from an ‘on’ and 4 from an ‘adjacent’. I also remove pins where camera malfunctions made this not possible (Pin 13 and Pin 2).

Removed unknowns, cryptic species, etc:

Since we could count some species as distinct individuals by life stages, I have combined those counts together here:

And built a little code to double check that it went correctly (should create two empty data frames)

I then deal with extreme schooling events following methods from the Donovan regimes paper: “Additional methodology was developed for dealing with outliers in the fish data, accounting for extreme observations of schooling species. Extreme observations in the database were defined by calculating the upper 99.9% of all individual observations (e.g. one species, size and count on an individual transect), resulting in 26 observations out of over 0.5 million, comprised of 11 species. The distribution of individual counts in the entire database for those 11 species was then used to identify observations that fell above the 99.0% quantile of counts for each species individually. These observations were adjusted to the 99.0% quantile for analysis.”

I want to make sure blocks are summarized appropriately, so I renamed them:

and ending by renaming the dataframe as just fish_tidy for future use: this allows me to make extra tidy data adjustments without breaking the dataframe name later on in the code!

Resource Fish Only

I’m pulling a list of fished genus from the paper “Perceptions and responses of Pacific Island fishers to changing coral reefs” by Rassweiler et al., 2021. I am editing fish_tidy to split genus and species, creating a column based off whether the genus is on the list for Rassweiler et al, and then making that into a binary variable ‘Resource’.

Also in this paper, it shows that “Naso, Chlorurus and Scarus collectively composed the bulk of the fishable biomass on the reef (48–66%) and a roughly similar total proportion of the catch (43–65%).” I have made a variable ‘most_targeted’ for these genus.

Wide data for ordination:

I’ve also made a wide version of the data frame for ordination:

Then used that one to re-make the long data frame with a presence column instead of abundance

Data Vis

Linda from Johanssen Lab recommended to make a plot that is ordered by the most abundant fish near the seep and look for any patterns, so I went and made this plot a number of ways. First, I kept all the species and looked a just the average abundance at all pins, to see if the communities looked relatively similar in this way.

All species average abundance

with Facet:

All species average abundance across gradient

This is too busy to read, and no major patterns jump out. I subset this to just the top 20 species along the gradient, and again the order of the pins does not have biological relevance.

Most common 20 species average abundance across gradient

So, now I will order the pins by those most common at pin 14

This isn’t showing the most obvious patterns, but there were only a few species of fish that swam in the sand near the seep. So, I’ll reorder the pins by the distance to the seep

Most common species average abundance at each pin by distance to seep:

To see if maybe a nutrient parameter would be good descriptor, I do this also by the low tide mean silicate

Most common species average abundance across the gradient of low tide mean silicate

Most common species average abundance across the gradient of nutrient delivery

Nutrient delivery quantified as a pc axis I got from the low tide mean of major nutrients across the gradient.This axis captures 72.87298% of the variation in mean low tide TA, Phosphate, Silicate, and N + N. Higer numbers on this axis indicate higher values of these variables.

Alpha Diversity:

Calculating richness against gradient metrics:

Calculating diversity against gradient metrics:

NMDS Non-Dimensional Ordination

Code from Linda:

I’ve been going back and forth on the transformation here, and am currently going with hellinger instead of vegan’s default, the wisconsin. This is because, from my limited understanding, the hellinger better handles cases where species are not present in multiple sites (known as double zeros).

UPDATE, after meeting with Megan Beacuse this is really one site, and we’re thinking about fish movement patterns, we’re doing a species maximum standardization.

## Run 0 stress 0.2836939 
## Run 1 stress 0.2916344 
## Run 2 stress 0.2845832 
## Run 3 stress 0.2972572 
## Run 4 stress 0.2848333 
## Run 5 stress 0.2910121 
## Run 6 stress 0.2925329 
## Run 7 stress 0.2848494 
## Run 8 stress 0.2886083 
## Run 9 stress 0.2865459 
## Run 10 stress 0.292077 
## Run 11 stress 0.2972778 
## Run 12 stress 0.2842212 
## Run 13 stress 0.2859633 
## Run 14 stress 0.2837173 
## ... Procrustes: rmse 0.005263253  max resid 0.03378612 
## Run 15 stress 0.284955 
## Run 16 stress 0.2868522 
## Run 17 stress 0.2964804 
## Run 18 stress 0.2985795 
## Run 19 stress 0.2875361 
## Run 20 stress 0.2907227 
## Run 21 stress 0.2931239 
## Run 22 stress 0.2836477 
## ... New best solution
## ... Procrustes: rmse 0.004414033  max resid 0.03557407 
## Run 23 stress 0.291555 
## Run 24 stress 0.2849 
## Run 25 stress 0.2934875 
## Run 26 stress 0.3016706 
## Run 27 stress 0.2852722 
## Run 28 stress 0.2900425 
## Run 29 stress 0.284922 
## Run 30 stress 0.300006 
## Run 31 stress 0.2942664 
## Run 32 stress 0.2970538 
## Run 33 stress 0.2874389 
## Run 34 stress 0.3001613 
## Run 35 stress 0.2861068 
## Run 36 stress 0.2870868 
## Run 37 stress 0.2932525 
## Run 38 stress 0.2855885 
## Run 39 stress 0.2858133 
## Run 40 stress 0.2875956 
## Run 41 stress 0.2909685 
## Run 42 stress 0.2900801 
## Run 43 stress 0.2849454 
## Run 44 stress 0.2872858 
## Run 45 stress 0.2967271 
## Run 46 stress 0.3050606 
## Run 47 stress 0.2879497 
## Run 48 stress 0.2927707 
## Run 49 stress 0.2904693 
## Run 50 stress 0.2875084 
## Run 51 stress 0.2969556 
## Run 52 stress 0.289401 
## Run 53 stress 0.28732 
## Run 54 stress 0.2839506 
## ... Procrustes: rmse 0.01084617  max resid 0.08258761 
## Run 55 stress 0.2911906 
## Run 56 stress 0.2952159 
## Run 57 stress 0.291694 
## Run 58 stress 0.2864163 
## Run 59 stress 0.2911648 
## Run 60 stress 0.2974377 
## Run 61 stress 0.2861896 
## Run 62 stress 0.2859992 
## Run 63 stress 0.2971073 
## Run 64 stress 0.2875706 
## Run 65 stress 0.2856663 
## Run 66 stress 0.2844447 
## Run 67 stress 0.2837142 
## ... Procrustes: rmse 0.004031453  max resid 0.02926313 
## Run 68 stress 0.2939105 
## Run 69 stress 0.2909872 
## Run 70 stress 0.2877214 
## Run 71 stress 0.2860971 
## Run 72 stress 0.2854211 
## Run 73 stress 0.2862497 
## Run 74 stress 0.2923067 
## Run 75 stress 0.2874524 
## Run 76 stress 0.2851625 
## Run 77 stress 0.2838468 
## ... Procrustes: rmse 0.006062516  max resid 0.03615627 
## Run 78 stress 0.287476 
## Run 79 stress 0.2894317 
## Run 80 stress 0.2933673 
## Run 81 stress 0.2883369 
## Run 82 stress 0.2927796 
## Run 83 stress 0.2896126 
## Run 84 stress 0.2877203 
## Run 85 stress 0.2987992 
## Run 86 stress 0.285382 
## Run 87 stress 0.284883 
## Run 88 stress 0.2856067 
## Run 89 stress 0.2859064 
## Run 90 stress 0.2944096 
## Run 91 stress 0.2848255 
## Run 92 stress 0.2873571 
## Run 93 stress 0.2842416 
## Run 94 stress 0.2873912 
## Run 95 stress 0.2845958 
## Run 96 stress 0.2918279 
## Run 97 stress 0.2836645 
## ... Procrustes: rmse 0.001609693  max resid 0.01111718 
## Run 98 stress 0.3004752 
## Run 99 stress 0.2998456 
## Run 100 stress 0.28651 
## Run 101 stress 0.2844491 
## Run 102 stress 0.2837545 
## ... Procrustes: rmse 0.006087905  max resid 0.03595345 
## Run 103 stress 0.2849559 
## Run 104 stress 0.2855576 
## Run 105 stress 0.2865474 
## Run 106 stress 0.2876973 
## Run 107 stress 0.2842487 
## Run 108 stress 0.2910615 
## Run 109 stress 0.285358 
## Run 110 stress 0.2937688 
## Run 111 stress 0.2937625 
## Run 112 stress 0.2982079 
## Run 113 stress 0.2859542 
## Run 114 stress 0.2881038 
## Run 115 stress 0.2842019 
## Run 116 stress 0.288367 
## Run 117 stress 0.2970652 
## Run 118 stress 0.3037476 
## Run 119 stress 0.283657 
## ... Procrustes: rmse 0.00244671  max resid 0.02081873 
## Run 120 stress 0.297279 
## Run 121 stress 0.2861792 
## Run 122 stress 0.2842019 
## Run 123 stress 0.2836512 
## ... Procrustes: rmse 0.0008324899  max resid 0.007460458 
## ... Similar to previous best
## *** Best solution repeated 1 times

## 
## Call:
## metaMDS(comm = mat.dis_fish, distance = "bray", k = 2, trymax = 500) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     mat.dis_fish 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.2836477 
## Stress type 1, weak ties
## Best solution was repeated 1 time in 123 tries
## The best solution was from try 22 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: scores missing

##        NMDS1       NMDS2 Site
## 1 -0.6647477 -0.15521066    B
## 2 -0.6915840 -0.25473280    B
## 3 -0.5496115 -0.03903441    B
## 4 -0.6823871 -0.29875479    B
## 5 -0.7267234 -0.20617498    B
## 6 -0.5718438 -0.24005037    B
## $vectors
##                                   NMDS1    NMDS2     r2 Pr(>r)    
## Acanthurus triostegus           0.99987 -0.01600 0.0269  0.161    
## Caranx melampygus               0.89993  0.43604 0.0013  0.918    
## Chaetodon vagabundus           -0.58799 -0.80887 0.3127  0.001 ***
## Cheilinus chlorourus           -0.96948 -0.24517 0.1350  0.002 ** 
## Halichoeres trimaculatus       -0.99105  0.13351 0.3401  0.001 ***
## Mulloidichthys vanicolensis    -0.87830  0.47810 0.0335  0.090 .  
## Rhinecanthus aculeatus         -0.98879  0.14932 0.2708  0.001 ***
## Scarus psittacus               -0.94732  0.32028 0.1914  0.001 ***
## Stegastes fasciolatus          -1.00000 -0.00124 0.0575  0.007 ** 
## Stethojulis bandanensis        -1.00000 -0.00241 0.1428  0.001 ***
## Balistapus undulatus            0.33816  0.94109 0.2699  0.001 ***
## Mulloidichthys flavolineatus   -0.65237  0.75790 0.0265  0.169    
## Stegastes sp                    0.83751 -0.54642 0.5961  0.001 ***
## Thalassoma hardwicke            0.69622  0.71783 0.0719  0.006 ** 
## Cheilodipterus quinquelineatus -0.90449 -0.42650 0.1078  0.002 ** 
## Parupeneus barberinus          -0.95335 -0.30186 0.1137  0.001 ***
## Pomacentrus pavo               -0.86474 -0.50223 0.0695  0.007 ** 
## Ostorhinchus nigrofasciatus    -0.46579 -0.88490 0.0689  0.007 ** 
## Scarus ghobban                 -0.94868 -0.31624 0.0393  0.052 .  
## Cheilinus trilobatus            0.92138  0.38866 0.0995  0.002 ** 
## Ostracion meleagris             0.48832 -0.87266 0.0002  0.989    
## Siganus spinus                 -0.91314  0.40764 0.0953  0.002 ** 
## Crenimugil crenilabrus         -0.72715 -0.68648 0.0530  0.027 *  
## Pateobatis fai                 -0.26739 -0.96359 0.0011  0.934    
## Chlorurus spilurus              0.20734  0.97827 0.2387  0.001 ***
## Echidna nebulosa               -0.74288 -0.66942 0.0271  0.163    
## Naso annulatus                 -0.61285 -0.79020 0.0272  0.161    
## Siganus argenteus              -0.61285 -0.79020 0.0272  0.161    
## Acanthurus nigrofuscus         -0.07154  0.99744 0.3601  0.001 ***
## Chaetodon ephippium             0.96173 -0.27399 0.0428  0.052 .  
## Cheilio inermis                -0.69530  0.71872 0.0543  0.022 *  
## Chromis viridis                -0.06684  0.99776 0.0429  0.048 *  
## Coris aygula                   -0.01420  0.99990 0.0080  0.605    
## Gomphosus varius                0.98002 -0.19892 0.2130  0.001 ***
## Parrotfish Unknonwn            -0.10411  0.99457 0.0086  0.570    
## Abudefduf sexfasciatus          0.57350  0.81921 0.0614  0.016 *  
## Parupeneus multifasciatus       0.63513  0.77241 0.0308  0.140    
## Aulostomus chinensis           -0.58804  0.80883 0.0051  0.705    
## Epinephelus merra               0.77469  0.63234 0.0080  0.583    
## Centropyge flavissima           0.38587  0.92255 0.1560  0.001 ***
## Oxycheilinus unifasciatus      -0.03342  0.99944 0.0276  0.170    
## Stegastes albifasciatus         0.28359  0.95894 0.0607  0.021 *  
## Zanclus cornutus                0.33056  0.94378 0.0559  0.019 *  
## Zebrasoma scopas                0.71106  0.70313 0.1346  0.001 ***
## Chaetodon citrinellus           0.55663  0.83076 0.0843  0.004 ** 
## Parupeneus insularis           -0.20349  0.97908 0.0162  0.352    
## Epibulus insidiator             0.99943 -0.03382 0.0930  0.001 ***
## Chaetodon auriga                0.97942  0.20183 0.0588  0.021 *  
## Naso lituratus                 -0.46236  0.88669 0.0156  0.367    
## Naso unicornis                  0.19704  0.98040 0.0065  0.738    
## Lethrinus olivaceus             0.49855  0.86686 0.0375  0.095 .  
## Chaetodon ulietensis            0.99977 -0.02161 0.0251  0.191    
## Dascyllus aruanus               0.29998 -0.95395 0.0699  0.005 ** 
## Carcharhinus melanopterus       0.56398  0.82579 0.0300  0.128    
## Chaetodon trifascialis         -0.30795  0.95140 0.0210  0.248    
## Chaetodon lunulatus             0.99113  0.13291 0.0435  0.046 *  
## Labroides dimidiatus            0.83435 -0.55124 0.1991  0.001 ***
## Abudefduf septemfasciatus       0.30670 -0.95181 0.0180  0.314    
## Chaetodon lunula                0.99877  0.04950 0.1384  0.001 ***
## Arothron meleagris              0.95754 -0.28829 0.0495  0.022 *  
## Hipposcarus longiceps           0.91383 -0.40610 0.0293  0.120    
## Parupeneus ciliatus            -0.30153  0.95346 0.0249  0.166    
## Aluterus scriptus               0.46712 -0.88420 0.0042  0.756    
## Cantherhines dumerilii          0.70586 -0.70835 0.0849  0.001 ***
## Scarus forsteni                 0.91300  0.40796 0.0047  0.753    
## Chrysiptera brownriggii         0.89177  0.45249 0.0485  0.038 *  
## Dascyllus trimaculatus          0.59216  0.80582 0.0221  0.232    
## Labroides bicolor               0.94925  0.31453 0.0146  0.416    
## Acanthurus guttatus             0.46701  0.88425 0.0073  0.652    
## Ostracion cubicus              -0.89469  0.44669 0.0030  0.833    
## Lutjanus fulvus                 0.05565  0.99845 0.0006  0.981    
## Abudefduf sordidus             -0.04007  0.99920 0.0270  0.191    
## Bodianus perditio              -0.36966  0.92917 0.0158  0.357    
## Chaetodon ornatissimus          0.20059 -0.97967 0.0074  0.626    
## Ellochelon vaigiensis           0.60308  0.79768 0.0038  0.801    
## Pseudocheilinus evanidus        0.40169  0.91578 0.0053  0.715    
## Zebrasoma velifer               0.73377  0.67940 0.0060  0.730    
## Pseudocheilinus tetrataenia     0.61457  0.78886 0.0021  0.860    
## Melichthys niger                0.40712  0.91337 0.0260  0.195    
## Fistularia commersonii          0.50824  0.86122 0.0172  0.333    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## $factors
## NULL
## 
## $na.action
## function (object, ...) 
## UseMethod("na.action")
## <bytecode: 0x114016b00>
## <environment: namespace:stats>
##                               NMDS1         NMDS2                  Species
## Chaetodon vagabundus     -0.3288165 -0.4523354776     Chaetodon vagabundus
## Cheilinus chlorourus     -0.3562309 -0.0900883972     Cheilinus chlorourus
## Halichoeres trimaculatus -0.5779498  0.0778611824 Halichoeres trimaculatus
## Rhinecanthus aculeatus   -0.5145087  0.0776953621   Rhinecanthus aculeatus
## Scarus psittacus         -0.4144708  0.1401281856         Scarus psittacus
## Stegastes fasciolatus    -0.2397427 -0.0002968729    Stegastes fasciolatus
##                           pval
## Chaetodon vagabundus     0.001
## Cheilinus chlorourus     0.002
## Halichoeres trimaculatus 0.001
## Rhinecanthus aculeatus   0.001
## Scarus psittacus         0.001
## Stegastes fasciolatus    0.007

##                                      NMDS1         NMDS2
## Chaetodon vagabundus           -0.32881651 -0.4523354776
## Cheilinus chlorourus           -0.35623094 -0.0900883972
## Halichoeres trimaculatus       -0.57794982  0.0778611824
## Rhinecanthus aculeatus         -0.51450873  0.0776953621
## Scarus psittacus               -0.41447084  0.1401281856
## Stegastes fasciolatus          -0.23974272 -0.0002968729
## Stethojulis bandanensis        -0.37789423 -0.0009120762
## Balistapus undulatus            0.17567443  0.4889035940
## Stegastes sp                    0.64663295 -0.4218891195
## Thalassoma hardwicke            0.18667744  0.1924696014
## Cheilodipterus quinquelineatus -0.29700786 -0.1400487619
## Parupeneus barberinus          -0.32147348 -0.1017893420
## Pomacentrus pavo               -0.22799546 -0.1324165965
## Ostorhinchus nigrofasciatus    -0.12225102 -0.2322506853
## Cheilinus trilobatus            0.29056538  0.1225671781
## Siganus spinus                 -0.28182491  0.1258107230
## Chlorurus spilurus              0.10130311  0.4779617250
## Acanthurus nigrofuscus         -0.04293312  0.5985638609
## Gomphosus varius                0.45233329 -0.0918122552
## Centropyge flavissima           0.15240946  0.3643893202
## Zebrasoma scopas                0.26088785  0.2579800960
## Chaetodon citrinellus           0.16160595  0.2411955237
## Epibulus insidiator             0.30478059 -0.0103129842
## Dascyllus aruanus               0.07930623 -0.2521997896
## Labroides dimidiatus            0.37229559 -0.2459715319
## Chaetodon lunula                0.37157904  0.0184170575
## Cantherhines dumerilii          0.20570134 -0.2064252692
##                                                       Species  pval
## Chaetodon vagabundus                     Chaetodon vagabundus 0.001
## Cheilinus chlorourus                     Cheilinus chlorourus 0.002
## Halichoeres trimaculatus             Halichoeres trimaculatus 0.001
## Rhinecanthus aculeatus                 Rhinecanthus aculeatus 0.001
## Scarus psittacus                             Scarus psittacus 0.001
## Stegastes fasciolatus                   Stegastes fasciolatus 0.007
## Stethojulis bandanensis               Stethojulis bandanensis 0.001
## Balistapus undulatus                     Balistapus undulatus 0.001
## Stegastes sp                                     Stegastes sp 0.001
## Thalassoma hardwicke                     Thalassoma hardwicke 0.006
## Cheilodipterus quinquelineatus Cheilodipterus quinquelineatus 0.002
## Parupeneus barberinus                   Parupeneus barberinus 0.001
## Pomacentrus pavo                             Pomacentrus pavo 0.007
## Ostorhinchus nigrofasciatus       Ostorhinchus nigrofasciatus 0.007
## Cheilinus trilobatus                     Cheilinus trilobatus 0.002
## Siganus spinus                                 Siganus spinus 0.002
## Chlorurus spilurus                         Chlorurus spilurus 0.001
## Acanthurus nigrofuscus                 Acanthurus nigrofuscus 0.001
## Gomphosus varius                             Gomphosus varius 0.001
## Centropyge flavissima                   Centropyge flavissima 0.001
## Zebrasoma scopas                             Zebrasoma scopas 0.001
## Chaetodon citrinellus                   Chaetodon citrinellus 0.004
## Epibulus insidiator                       Epibulus insidiator 0.001
## Dascyllus aruanus                           Dascyllus aruanus 0.005
## Labroides dimidiatus                     Labroides dimidiatus 0.001
## Chaetodon lunula                             Chaetodon lunula 0.001
## Cantherhines dumerilii                 Cantherhines dumerilii 0.001

Constrained Ordination - Distance based redundancy analysis

I started with using the pulse pc axis, the low tide mean silicate, the distance to seep, and the distance to shore:

Tests revealed the shore distance didn’t matter (permutation test p = 0.555), so I changed this to only use the pulse pc axis, the low tide mean silicate, the distance to seep, and the distance to shore:

## Call: capscale(formula = transformed_communities ~
## Low_Tide_Mean_Phosphate_umolL + Low_Tide_Mean_Silicate_umolL +
## dist_to_seep_m + Low_Tide_Mean_NN_umolL + Low_Tide_Mean_Ammonia_umolL,
## data = transformed_communities_with_explan, distance = "bray")
## 
##               Inertia Proportion Rank
## Total         26.8723                
## RealTotal     34.6230     1.0000     
## Constrained    4.4752     0.1293    5
## Unconstrained 30.1478     0.8707   55
## Imaginary     -7.7507                
## Inertia is squared Bray distance 
## Species scores projected from 'transformed_communities' 
## 
## Eigenvalues for constrained axes:
##   CAP1   CAP2   CAP3   CAP4   CAP5 
## 2.2677 1.0523 0.6100 0.3174 0.2278 
## 
## Eigenvalues for unconstrained axes:
##  MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8 
## 4.167 2.411 1.965 1.709 1.645 1.407 1.325 1.175 
## (Showing 8 of 55 unconstrained eigenvalues)

The output reports the total inertia, which is the total amount of variation (dissimilarity) in the data. This inertia is decomposed into ‘constrained’ and ‘unconstrained’ components. The constrained component is the total amount of variation explained by the predictors (18.36%), while the unconstrained component is the remaining ‘residual’ variation. There is also info on ‘real’ and ‘imaginary’ components, due to the negative eigenvalues issue which arises with PCoA.

The default plot shows how these variables are loaded onto the first two CAP axes, and shows how the samples (circles) are ordinated on those axes, as well as the species scores (red crosses).

We can get the variance explained by these axes from summary:

## $importance
## Importance of components:
##                         CAP1    CAP2    CAP3     CAP4     CAP5   MDS1    MDS2
## Eigenvalue            2.2677 1.05234 0.60998 0.317403 0.227793 4.1670 2.41098
## Proportion Explained  0.0655 0.03039 0.01762 0.009167 0.006579 0.1204 0.06964
## Cumulative Proportion 0.0655 0.09589 0.11351 0.122676 0.129255 0.2496 0.31924
##                          MDS3    MDS4    MDS5    MDS6    MDS7    MDS8   MDS9
## Eigenvalue            1.96497 1.70933 1.64517 1.40655 1.32534 1.17545 1.0524
## Proportion Explained  0.05675 0.04937 0.04752 0.04062 0.03828 0.03395 0.0304
## Cumulative Proportion 0.37600 0.42537 0.47288 0.51351 0.55179 0.58574 0.6161
##                         MDS10   MDS11   MDS12   MDS13   MDS14   MDS15   MDS16
## Eigenvalue            1.02505 0.94062 0.81301 0.78152 0.72257 0.71529 0.67685
## Proportion Explained  0.02961 0.02717 0.02348 0.02257 0.02087 0.02066 0.01955
## Cumulative Proportion 0.64574 0.67291 0.69639 0.71896 0.73983 0.76049 0.78004
##                         MDS17   MDS18   MDS19   MDS20   MDS21   MDS22   MDS23
## Eigenvalue            0.64167 0.60873 0.51663 0.47472 0.46950 0.42736 0.42173
## Proportion Explained  0.01853 0.01758 0.01492 0.01371 0.01356 0.01234 0.01218
## Cumulative Proportion 0.79857 0.81615 0.83107 0.84479 0.85835 0.87069 0.88287
##                         MDS24    MDS25    MDS26    MDS27    MDS28    MDS29
## Eigenvalue            0.37365 0.338154 0.332687 0.300045 0.254459 0.249022
## Proportion Explained  0.01079 0.009767 0.009609 0.008666 0.007349 0.007192
## Cumulative Proportion 0.89366 0.903428 0.913037 0.921703 0.929053 0.936245
##                          MDS30    MDS31    MDS32   MDS33  MDS34    MDS35
## Eigenvalue            0.229672 0.217394 0.198832 0.17348 0.1662 0.160803
## Proportion Explained  0.006634 0.006279 0.005743 0.00501 0.0048 0.004644
## Cumulative Proportion 0.942879 0.949157 0.954900 0.95991 0.9647 0.969355
##                          MDS36   MDS37    MDS38    MDS39    MDS40    MDS41
## Eigenvalue            0.137132 0.11600 0.108065 0.094114 0.085650 0.079762
## Proportion Explained  0.003961 0.00335 0.003121 0.002718 0.002474 0.002304
## Cumulative Proportion 0.973315 0.97667 0.979787 0.982505 0.984979 0.987283
##                          MDS42    MDS43    MDS44    MDS45    MDS46    MDS47
## Eigenvalue            0.075081 0.066099 0.056860 0.046630 0.038583 0.035702
## Proportion Explained  0.002169 0.001909 0.001642 0.001347 0.001114 0.001031
## Cumulative Proportion 0.989451 0.991360 0.993003 0.994350 0.995464 0.996495
##                           MDS48     MDS49     MDS50     MDS51     MDS52
## Eigenvalue            0.0307919 0.0265469 0.0191956 0.0148918 0.0107622
## Proportion Explained  0.0008893 0.0007667 0.0005544 0.0004301 0.0003108
## Cumulative Proportion 0.9973844 0.9981511 0.9987056 0.9991357 0.9994465
##                           MDS53     MDS54     MDS55
## Eigenvalue            0.0089720 0.0078800 2.311e-03
## Proportion Explained  0.0002591 0.0002276 6.676e-05
## Cumulative Proportion 0.9997056 0.9999332 1.000e+00

The first CAP axis explains 6.97% of the constrained community variation, and the second axis explains 4.34%. Therefore, this represents 2.07% of the total community variation.

Though this is already pretty poor, I might as well finish it out as a coding exercise. The loading coefficients for the explanatory variables on the constrained axes are:

##                                      CAP1        CAP2       CAP3        CAP4
## Low_Tide_Mean_Phosphate_umolL -0.14483502  0.38927986 -0.5113166 -0.34042879
## Low_Tide_Mean_Silicate_umolL   0.01084232  0.38145657 -0.8480870 -0.28435951
## dist_to_seep_m                 0.71378221 -0.60687018 -0.1822926  0.02825845
## Low_Tide_Mean_NN_umolL         0.29683934  0.05191737 -0.7186106 -0.62662481
## Low_Tide_Mean_Ammonia_umolL    0.52898679  0.62565553  0.0345422 -0.17859909
##                                     CAP5 MDS1
## Low_Tide_Mean_Phosphate_umolL -0.6709304    0
## Low_Tide_Mean_Silicate_umolL  -0.2329408    0
## dist_to_seep_m                -0.2969755    0
## Low_Tide_Mean_NN_umolL        -0.0114547    0
## Low_Tide_Mean_Ammonia_umolL   -0.5437254    0

To test the important of all the predictors in combination:

## Permutation test for capscale under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = transformed_communities ~ Low_Tide_Mean_Phosphate_umolL + Low_Tide_Mean_Silicate_umolL + dist_to_seep_m + Low_Tide_Mean_NN_umolL + Low_Tide_Mean_Ammonia_umolL, data = transformed_communities_with_explan, distance = "bray")
##           Df SumOfSqs      F Pr(>F)    
## Model      5   4.4752 3.8595  0.001 ***
## Residual 130  30.1478                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

So this tests whether the total variation explained by the constrained axes is significant. Interestingly, though the amount of variation seems trivial, it is showing up as significant…

## Permutation test for capscale under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = transformed_communities ~ Low_Tide_Mean_Phosphate_umolL + Low_Tide_Mean_Silicate_umolL + dist_to_seep_m + Low_Tide_Mean_NN_umolL + Low_Tide_Mean_Ammonia_umolL, data = transformed_communities_with_explan, distance = "bray")
##           Df SumOfSqs      F Pr(>F)    
## CAP1       1   2.2677 9.7785  0.001 ***
## CAP2       1   1.0523 4.5378  0.001 ***
## CAP3       1   0.6100 2.6303  0.005 ** 
## CAP4       1   0.3174 1.3687  0.283    
## CAP5       1   0.2278 0.9823  0.450    
## Residual 130  30.1478                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

So, both the axes we have plotted are considered significant.

## Permutation test for capscale under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = transformed_communities ~ Low_Tide_Mean_Phosphate_umolL + Low_Tide_Mean_Silicate_umolL + dist_to_seep_m + Low_Tide_Mean_NN_umolL + Low_Tide_Mean_Ammonia_umolL, data = transformed_communities_with_explan, distance = "bray")
##                                Df SumOfSqs      F Pr(>F)    
## Low_Tide_Mean_Phosphate_umolL   1   0.8001 3.4501  0.001 ***
## Low_Tide_Mean_Silicate_umolL    1   0.5079 2.1902  0.006 ** 
## dist_to_seep_m                  1   0.8770 3.7819  0.001 ***
## Low_Tide_Mean_NN_umolL          1   0.3461 1.4926  0.071 .  
## Low_Tide_Mean_Ammonia_umolL     1   1.5439 6.6574  0.001 ***
## Residual                      130  30.1478                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interestingly, all the axes are coming up as significant contributors to this variation.

PERMANOVA:

## Permutation test for adonis under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = transformed_communities ~ Low_Tide_Mean_Phosphate_umolL + Low_Tide_Mean_Silicate_umolL + dist_to_seep_m + Low_Tide_Mean_NN_umolL + Low_Tide_Mean_Ammonia_umolL, data = transformed_communities_with_explan, by = "margin", dist = "bray")
##                                Df SumOfSqs      R2      F Pr(>F)    
## Low_Tide_Mean_Phosphate_umolL   1   0.7521 0.02799 4.3231  0.001 ***
## Low_Tide_Mean_Silicate_umolL    1   0.4638 0.01726 2.6661  0.010 ** 
## dist_to_seep_m                  1   0.8373 0.03116 4.8128  0.001 ***
## Low_Tide_Mean_NN_umolL          1   0.2887 0.01074 1.6593  0.070 .  
## Low_Tide_Mean_Ammonia_umolL     1   1.5027 0.05592 8.6376  0.001 ***
## Residual                      130  22.6158 0.84160                  
## Total                         135  26.8723 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Again, all of our variables are significant contributors to community composition, though with very low r2.

Mixed Models for Multivariate Data

The first model I test is looking at the probability of presence, and whether it differs between the CV of Salinity. With this, the effect of salinity differs randomly around the species. The random effect for salinity can quantify how much the community composition chnanges along the gradient. If the random effect is non-0, the species change in relative abundance along the salinity gradient.

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Presence ~ CV_Salinity + (1 + CV_Salinity | species)
##    Data: presence_explan
## 
##      AIC      BIC   logLik deviance df.resid 
##   6442.4   6478.9  -3216.2   6432.4    10875 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4621 -0.3220 -0.1567 -0.1048  9.5805 
## 
## Random effects:
##  Groups  Name        Variance  Std.Dev. Corr 
##  species (Intercept)     5.338   2.31        
##          CV_Salinity 56749.425 238.22   -0.56
## Number of obs: 10880, groups:  species, 80
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.3271     0.2262  -10.29   <2e-16 ***
## CV_Salinity -86.2975     8.2852  -10.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## CV_Salinity -0.158
## $species

We can test the salinity effect by dropping the random slope and doing an LRT

## Data: presence_explan
## Models:
## mod.no.salinity: Presence ~ CV_Salinity + (1 | species)
## mod.salinity: Presence ~ CV_Salinity + (1 + CV_Salinity | species)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod.no.salinity    3 6478.2 6500.0 -3236.1   6472.2                         
## mod.salinity       5 6442.4 6478.9 -3216.2   6432.4 39.769  2  2.314e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This doesn’t look great, but the technique is interesting.

GAM

I want to try to do this to look at non-linear distribution along the gradient. I think some species might be non-linear as they might avoid the highly fresh region. I think this might be easiest to look at if I only looked at certain species of interest, which I still would like to define.

For now, I’m setting up the code using species driving differences in the nmds

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## Presence ~ s(dist_to_seep_m, by = species)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.33048    0.04101  -32.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                                           edf Ref.df Chi.sq
## s(dist_to_seep_m):speciesAbudefduf sexfasciatus         1.000  1.000  0.757
## s(dist_to_seep_m):speciesAcanthurus nigrofuscus         8.815  8.971 55.883
## s(dist_to_seep_m):speciesArothron meleagris             1.000  1.000  0.015
## s(dist_to_seep_m):speciesBalistapus undulatus           5.940  7.046 24.486
## s(dist_to_seep_m):speciesCantherhines dumerilii         1.000  1.000  0.001
## s(dist_to_seep_m):speciesCentropyge flavissima          8.459  8.904 30.798
## s(dist_to_seep_m):speciesChaetodon auriga               7.703  8.539 29.454
## s(dist_to_seep_m):speciesChaetodon citrinellus          6.397  7.484 43.181
## s(dist_to_seep_m):speciesChaetodon lunula               2.295  2.859 29.269
## s(dist_to_seep_m):speciesChaetodon lunulatus            1.160  1.302  1.720
## s(dist_to_seep_m):speciesChaetodon vagabundus           5.472  6.459 29.711
## s(dist_to_seep_m):speciesCheilinus chlorourus           7.811  8.559 36.293
## s(dist_to_seep_m):speciesCheilinus trilobatus           1.757  2.191 13.263
## s(dist_to_seep_m):speciesCheilio inermis                2.472  3.072  6.393
## s(dist_to_seep_m):speciesCheilodipterus quinquelineatus 1.449  1.764  1.745
## s(dist_to_seep_m):speciesChlorurus spilurus             8.306  8.850 58.604
## s(dist_to_seep_m):speciesChromis viridis                6.429  7.468 14.236
## s(dist_to_seep_m):speciesChrysiptera brownriggii        1.000  1.000  1.636
## s(dist_to_seep_m):speciesCrenimugil crenilabrus         1.000  1.000  0.447
## s(dist_to_seep_m):speciesDascyllus aruanus              1.000  1.000  0.887
## s(dist_to_seep_m):speciesEpibulus insidiator            1.000  1.000  1.573
## s(dist_to_seep_m):speciesGomphosus varius               5.972  7.075 41.872
## s(dist_to_seep_m):speciesHalichoeres trimaculatus       3.792  4.604 18.778
## s(dist_to_seep_m):speciesLabroides dimidiatus           1.949  2.427  1.688
## s(dist_to_seep_m):speciesOstorhinchus nigrofasciatus    1.000  1.000  0.709
## s(dist_to_seep_m):speciesParupeneus barberinus          1.579  1.944  4.884
## s(dist_to_seep_m):speciesPomacentrus pavo               1.000  1.000  0.455
## s(dist_to_seep_m):speciesRhinecanthus aculeatus         2.964  3.647  3.785
## s(dist_to_seep_m):speciesScarus psittacus               8.953  8.999 52.759
## s(dist_to_seep_m):speciesSiganus spinus                 1.801  2.237  4.023
## s(dist_to_seep_m):speciesStegastes albifasciatus        1.000  1.000  0.332
## s(dist_to_seep_m):speciesStegastes fasciolatus          1.000  1.000  0.776
## s(dist_to_seep_m):speciesStegastes sp                   1.264  1.480  1.514
## s(dist_to_seep_m):speciesStethojulis bandanensis        8.952  8.999 57.191
## s(dist_to_seep_m):speciesThalassoma hardwicke           6.545  7.610 84.057
## s(dist_to_seep_m):speciesZanclus cornutus               1.000  1.001  0.541
## s(dist_to_seep_m):speciesZebrasoma scopas               6.935  7.956 31.203
##                                                          p-value    
## s(dist_to_seep_m):speciesAbudefduf sexfasciatus         0.384341    
## s(dist_to_seep_m):speciesAcanthurus nigrofuscus          < 2e-16 ***
## s(dist_to_seep_m):speciesArothron meleagris             0.904621    
## s(dist_to_seep_m):speciesBalistapus undulatus           0.000964 ***
## s(dist_to_seep_m):speciesCantherhines dumerilii         0.979198    
## s(dist_to_seep_m):speciesCentropyge flavissima          0.000875 ***
## s(dist_to_seep_m):speciesChaetodon auriga               0.000428 ***
## s(dist_to_seep_m):speciesChaetodon citrinellus           < 2e-16 ***
## s(dist_to_seep_m):speciesChaetodon lunula               2.19e-06 ***
## s(dist_to_seep_m):speciesChaetodon lunulatus            0.311018    
## s(dist_to_seep_m):speciesChaetodon vagabundus           7.05e-05 ***
## s(dist_to_seep_m):speciesCheilinus chlorourus           3.33e-05 ***
## s(dist_to_seep_m):speciesCheilinus trilobatus           0.001728 ** 
## s(dist_to_seep_m):speciesCheilio inermis                0.103212    
## s(dist_to_seep_m):speciesCheilodipterus quinquelineatus 0.265255    
## s(dist_to_seep_m):speciesChlorurus spilurus              < 2e-16 ***
## s(dist_to_seep_m):speciesChromis viridis                0.056347 .  
## s(dist_to_seep_m):speciesChrysiptera brownriggii        0.200933    
## s(dist_to_seep_m):speciesCrenimugil crenilabrus         0.503831    
## s(dist_to_seep_m):speciesDascyllus aruanus              0.346236    
## s(dist_to_seep_m):speciesEpibulus insidiator            0.209870    
## s(dist_to_seep_m):speciesGomphosus varius                < 2e-16 ***
## s(dist_to_seep_m):speciesHalichoeres trimaculatus       0.001729 ** 
## s(dist_to_seep_m):speciesLabroides dimidiatus           0.438186    
## s(dist_to_seep_m):speciesOstorhinchus nigrofasciatus    0.399672    
## s(dist_to_seep_m):speciesParupeneus barberinus          0.060421 .  
## s(dist_to_seep_m):speciesPomacentrus pavo               0.499865    
## s(dist_to_seep_m):speciesRhinecanthus aculeatus         0.363621    
## s(dist_to_seep_m):speciesScarus psittacus                < 2e-16 ***
## s(dist_to_seep_m):speciesSiganus spinus                 0.143617    
## s(dist_to_seep_m):speciesStegastes albifasciatus        0.564684    
## s(dist_to_seep_m):speciesStegastes fasciolatus          0.378432    
## s(dist_to_seep_m):speciesStegastes sp                   0.443667    
## s(dist_to_seep_m):speciesStethojulis bandanensis         < 2e-16 ***
## s(dist_to_seep_m):speciesThalassoma hardwicke            < 2e-16 ***
## s(dist_to_seep_m):speciesZanclus cornutus               0.462383    
## s(dist_to_seep_m):speciesZebrasoma scopas               0.000126 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.197   Deviance explained = 13.5%
## UBRE = 0.059816  Scale est. = 1         n = 5032

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## Presence ~ s(dist_to_seep_m, by = species)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.9396     0.0641  -30.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                                         edf Ref.df  Chi.sq
## s(dist_to_seep_m):speciesAcanthurus guttatus          1.000  1.000   0.050
## s(dist_to_seep_m):speciesAcanthurus nigrofuscus       8.813  8.971  90.316
## s(dist_to_seep_m):speciesAcanthurus triostegus        8.940  8.999  70.134
## s(dist_to_seep_m):speciesCaranx melampygus            1.000  1.000   0.000
## s(dist_to_seep_m):speciesCheilinus chlorourus         8.015  8.635  55.889
## s(dist_to_seep_m):speciesCheilinus trilobatus         8.848  8.989  46.121
## s(dist_to_seep_m):speciesChlorurus spilurus           8.438  8.900 104.930
## s(dist_to_seep_m):speciesEpinephelus merra            2.999  3.740   4.114
## s(dist_to_seep_m):speciesLethrinus olivaceus          2.838  3.535   5.633
## s(dist_to_seep_m):speciesLutjanus fulvus              1.000  1.000   0.008
## s(dist_to_seep_m):speciesMulloidichthys flavolineatus 4.704  5.722  18.225
## s(dist_to_seep_m):speciesMulloidichthys vanicolensis  1.752  2.176   1.530
## s(dist_to_seep_m):speciesNaso annulatus               1.000  1.000   0.168
## s(dist_to_seep_m):speciesNaso lituratus               2.629  3.250   1.347
## s(dist_to_seep_m):speciesNaso unicornis               1.000  1.000   0.057
## s(dist_to_seep_m):speciesParupeneus barberinus        1.619  1.996   8.281
## s(dist_to_seep_m):speciesParupeneus ciliatus          1.000  1.000   0.003
## s(dist_to_seep_m):speciesParupeneus insularis         1.000  1.000   0.437
## s(dist_to_seep_m):speciesParupeneus multifasciatus    7.643  8.485  76.935
## s(dist_to_seep_m):speciesScarus forsteni              1.000  1.000   0.148
## s(dist_to_seep_m):speciesScarus ghobban               1.232  1.426   2.025
## s(dist_to_seep_m):speciesScarus psittacus             8.993  9.000  92.985
## s(dist_to_seep_m):speciesSiganus argenteus            1.000  1.000   0.168
## s(dist_to_seep_m):speciesSiganus spinus               8.165  8.675  25.875
##                                                       p-value    
## s(dist_to_seep_m):speciesAcanthurus guttatus          0.82283    
## s(dist_to_seep_m):speciesAcanthurus nigrofuscus       < 2e-16 ***
## s(dist_to_seep_m):speciesAcanthurus triostegus        < 2e-16 ***
## s(dist_to_seep_m):speciesCaranx melampygus            0.99777    
## s(dist_to_seep_m):speciesCheilinus chlorourus         < 2e-16 ***
## s(dist_to_seep_m):speciesCheilinus trilobatus         < 2e-16 ***
## s(dist_to_seep_m):speciesChlorurus spilurus           < 2e-16 ***
## s(dist_to_seep_m):speciesEpinephelus merra            0.33378    
## s(dist_to_seep_m):speciesLethrinus olivaceus          0.16821    
## s(dist_to_seep_m):speciesLutjanus fulvus              0.92724    
## s(dist_to_seep_m):speciesMulloidichthys flavolineatus 0.00473 ** 
## s(dist_to_seep_m):speciesMulloidichthys vanicolensis  0.54852    
## s(dist_to_seep_m):speciesNaso annulatus               0.68157    
## s(dist_to_seep_m):speciesNaso lituratus               0.68487    
## s(dist_to_seep_m):speciesNaso unicornis               0.81216    
## s(dist_to_seep_m):speciesParupeneus barberinus        0.01382 *  
## s(dist_to_seep_m):speciesParupeneus ciliatus          0.95788    
## s(dist_to_seep_m):speciesParupeneus insularis         0.50870    
## s(dist_to_seep_m):speciesParupeneus multifasciatus    < 2e-16 ***
## s(dist_to_seep_m):speciesScarus forsteni              0.70021    
## s(dist_to_seep_m):speciesScarus ghobban               0.17823    
## s(dist_to_seep_m):speciesScarus psittacus             < 2e-16 ***
## s(dist_to_seep_m):speciesSiganus argenteus            0.68157    
## s(dist_to_seep_m):speciesSiganus spinus               0.00329 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.261   Deviance explained = 18.1%
## UBRE = -0.12018  Scale est. = 1         n = 3264

This code didn’t really work as intended. But, when I talked to Kyle Edwards, he suggested that I make a different gam for each species. This is more or less a preliminary step to see if a linear fit is good enough for this data:

From further inspection, the model has issues with Naso annulatus, Naso unicornis, Scarus oviceps

print(warnings_list[[“Naso unicornis”]]) [1] “Iteration limit reached without full convergence - check carefully”

print(warnings_list[[“Naso annulatus”]]) [1] “Fitting terminated with step failure - check results carefully”

print(warnings_list[[“Scarus oviceps”]]) [1] “Fitting terminated with step failure - check results carefully”